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Abstract 

The influence of viscosity contrast on buoyantly unstable miscible fluids in a porous medium is 
investigated through a linear stability analysis (LSA) as well as direct numerical simulations (DNS). 
The linear stability method implemented in this paper is based on an initial value approach, which 
helps to capture the onset of instability more accurately than the quasi-steady state analysis. In 
the absence of displacement, we show that viscosity contrast delays the onset of instability in 
buoyantly unstable miscible fluids. Further, it is observed that suitably choosing the viscosity 
contrast and injection velocity a gravitationally unstable miscible interface can be stabilised com¬ 
pletely. Through LSA we draw a phase diagram, which shows three distinct stability regions in a 
parameter space spanned by the displacement velocity and the viscosity contrast. DNS are per¬ 
formed corresponding to parameters from each regime and the results obtained are in accordance 
with the linear stability results. Moreover, the conversion from a dimensionless formulation to the 
other and its essence to compare between two different type of flow problems associated with each 
dimensionless formulation are discussed. 
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I. INTRODUCTION 


Understanding hydrodynamic instabilities and mixing of miscible fluids in porous media 
is an active area of research, related to several industrial and environmental processes, such 
as oil recovery |8], CO 2 sequestration m, groundwater contamination [2], chromatography 
[7], to name a few. Global warming attributed to anthropogenic emission of greenhouse 
gases is one of the major challenges facing mankind. Carbon capture and storage (or CO 2 
sequestration) in underground aquifers is detected as a promising mean to restrict unwanted 
rise of greenhouse gas level in the atmosphere. Geological sequestration of disposal in deep 
saline aquifers offers several key assets of miscible hngering instabilities driven by both 
viscous and buoyancy forces. At a subsurface saline aquifer site the injected super-critical 
lighter CO 2 rises up and accumulates under an impervious rock, followed by dissolution of 
CO 2 into the underlying brine. An unstably stratihed diffusive interface of heavier CO 2 
dissolved brine is formed above the pure brine and transition to natural convection in the 
form of unstable sinking plumes is featured in time. 

Several theoretical [lEllS] and experimental [DOS] studies are conducted to understand 
the convective instabilities devoted to characterising optimal storage sites that ensure CO 2 
does not leak into the environment. In such convective flows, viscosity variation, albeit small, 
at the diffusive interface influences the onset of instability and hence plays significant role in 
characterising the storage sites. Despite having enormous importance in real life applications 
this problem remains poorly explored. Manickam and Homsy [13] have shown through LSA 
as well as DNS that locally stable regions can be introduced by suitably choosing viscosity 
prohle and injection velocity to buoyantly unstable diffusive interface. Recently, Daniel 
and Riaz |3] used fixed interface and moving interface models to compare the theoretical 
predictions with the corresponding experimental observations [1]. These authors showed, 
through an LSA, that in the absence of displacement the onset time increases when the 
dynamic viscosity increases with the depth. On the other hand, when the more viscous 
fluid displaces the less viscous one from above the onset time depends non-monotonically 
on the viscosity contrast between the two fluids. In both these studies LSA was performed 
under a quasi-steady state approximation, which has its own drawback, since it does not 
capture possible transient growth of the linearly unstable modes and hence fails to predict 
the onset of instability accurately. For flows with unsteady base-state, transient growth of 
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deldx = 0, u = (U, 0 ) 



FIG. 1: Schematic of the flow. 


the perturbations is possible in hydrodynamic stability problems driven solely by buoyancy 
force [IH] or viscous force [9]. 

Discussion of the above-mentioned literature reveals that transient growth plays impor¬ 
tant role on the onset of instability. Therefore, it is essential to discuss an LSA without 
quasi-steady state approximation, which not only captures the onset time more accurately, 
but also represents the physics appropriately. Recently, Hota et al. [10] have discussed 
an LSA based on an initial value problem (IVP) approach using a Fourier pseudo-spectral 
method. This IVP based LSA method captures the diffusion dominated region, which was 
never captured before using quasi-steady state approximation. In this context, we present 
an LSA [10] and DNS [20] using a Fourier pseudo-spectral method to analyse the influence 
of viscosity contrast on buoyantly unstable miscible fluids in vertical porous media when 
the dynamic viscosity of the upper fluid is more or less than that of the lower fluid. From 
linear stability results it is identihed that in the absence of displacement, viscosity contrast 
of either kind delays the onset of instability of gravitationally unstable miscible fluids. We 
further identify two different dynamical regimes in which the instability is dominated either 
by the viscous force or the buoyancy force. Our results pave the way to new insights into 
the influence of viscosity contrast on hngering instability in a buoyantly unstable miscible 
system. 
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II. PROBLEM FORMULATION 


A. Governing equations 

Let us consider the displacement of a fluid of dynamic viscosity /iiow and density piow 
by another fluid of dynamic viscosity Pupp and density Pupp (> Plow) in a vertical porous 
medium, where the subscripts ‘low’ and ‘upp’ correspond to the lower and upper fluids, 
respectively. The fluids are assumed to be incompressible and miscible with each other. We 
further assume that the displacing fluid is injected at a uniform speed U vertically downward 
as shown in Fig. [T] The dynamic viscosity of the fluid varies with a solute concentration, 
c, which satisfies a convection-diffusion equation. The fluid velocity can be determined in 
terms of the Darcy’s law, which is a mathematical analogous to the flow equation in 2D 
homogeneous porous media. Hence, the governing equations can be written as. 


<1 

IS 

o 

(1) 

Vp- ^^K + p{c)g, 

(2) 

K - 

dtc + u- Vc = DV^c, 

(3) 


where k is the permeability of the porous media, g = {g, 0) with g being the gravitational 
acceleration and D is the isotropic dispersion co-efficient. 


B. Dimensionless formulation 

In order to obtain dimensionless equations we use Lch = Vch/D and Teh = as 

the characteristic length and time scales, respectively, where Wh = |Ap|Kp//ich is the buoy¬ 
ancy induced velocity, and pch is the characteristic viscosity. The characteristic pressure, 
concentration and density are taken to be Pch-D/^, Cupp and |Ap| = |pupp — piow|, respec¬ 
tively. The related dimensionless equations in a Lagrangian frame of reference moving with 
a dimensionless velocity, U = U/Vch, are 

V ■ M = 0 

Vp = -p(c) {u + Ue^) -F p(c)e„, 
dtc + u- Vc = V^c, 
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(4) 

(5) 

( 6 ) 



where ex is the unit vector in the x-direction. We assume that the density varies linearly 
with concentration, such that the dimensionless density prohle is given as p{c) = c, 

lAc) = ( 7 ) 


where /(c) is a linear function of c and the log-mobility ratio R represents the natural 
logarithm of the ratio of the dynamic viscosities of two fluids. The explicit form of /(c) and 
R depend on the characteristic viscosity /ich- We define /ich, and the corresponding /(c) and 


R appropriately in (IV while discussing the numerical results obtained from LSA and DNS. 


C. Initial and boundary conditions 


A description of appropriate initial and boundary conditions makes the mathematical 
formulation of the above problem complete. In the Lagrangian frame of reference the initial 
conditions for the concentration and velocity can be written as, 


c = 


1, for a: < 0 ^ , 

, and u = (0,0). 

0, for X > 0 


( 8 ) 


Along the longitudinal boundaries we have u = (0, 0) and dc/dx —)■ 0 as x —>■ ±cxd, while the 
transverse boundary conditions are dv/dy = Q (constant pressure cf. [IB]) and dc/dy —)■ 0. 


III. LINEAR STABILITY ANALYSIS 

In this section we discuss linear stability analysis of hngering instabilities driven by both 
viscosity and density contrasts in miscible displacements. 


A. Stream function formulation 

For a two dimensional flow the continuity equation can be satished identically by intro¬ 
ducing a stream function, '0(x,|/,t), such that u = d'lfj/dy and v = —dtlj/dx. Taking curl of 
equation ([^ and representing velocities in terms of the stream function, equations (|^ and 
0 are recast as, 

= -RVf{c) (Vc • + Udyc) + -^dyc, 

jX\C) 

dtc + {dy'ilj){dxc) - {dx'ilj){dyc) = V^c, 
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(9) 

( 10 ) 




where V = d/dc. 


B. Linearized perturbation equations 


For base-state flow we assume = 0 that implies "06 to be constant, which can be 
assumed to be equal to 0 without any loss of generality. We also assume that the base-state 
concentration is homogeneous in the ^/-direction, i.e. Cb = Cb{x,t). Under these assumptions 
base-state flow is given by decaying error function solution of step-like initial concentration 
profile, i.e., Cb{x,t) = 0.5 eric{x/2\/t). Introduce infinitesimal perturbations such that 
c{x,y,t) = Cb{x,t) + c'{x,y,t) and 'ip{x,y,t) = 'ipb + 'ip'{x,y,t), etc., and substitute these in 
equations ([^ and (10) to obtain. 


dtCb + dtc' + {d^Cb + dxd) dy'ip' - {d^'ilj'){dyc!) = V Cb + V c'. 


= -RVficb + c') [(Vc' + {d,Cb)ey) ■ + Udyc'] + 


—d d 

r.l\ ’ 


( 11 ) 

( 12 ) 


p(cfe d) 

with ey being the unit vector in the ^/-direction. Subtracting the base-state equations from 
equations ( [IT| ) and (12) we obtain the following coupled nonlinear partial differential equa¬ 
tions in terms of the perturbation quantities d and iIj' as, 

^ic' {d^Cb + d^d) dyd' - {dyd){d^'tp') = V^c', (13) 

V"V’' = -RVf{c) [(Vc' + {d,Cb)ey) ■ Wd' + Udyd] + (14) 

Linearizing these equations in terms of the perturbation quantities we obtain, 
did + dxCbdyij' = V^c', 


= -RVfic) [id^Cb)id,d’) + Udyd] + 


- d d 

’ ■ 


( 15 ) 

( 16 ) 


We solve equations (15) and (16) using a pseudo-spectral method to obtain the spatio- 


temporal evolution of the perturbation quantities, d{x, y, t) and V, t) and calculate the 
growth rates associated with concentration and velocity perturbations liniiig, 


1 dEr' 


1 dE„i 


(Jr. = 


2Ec’ dt ’ ^ 2Eyi dt ’ 

from the amplification measures defined as |10j . 

{dfdxdy, Ey> = [ [ {dyd'f + id^d'f dxdy, E = Er' + E^>. 


a = 


1 dE 


(17) 


Erl — 



Following Hota et al. [TO] we have used a (defined in equation (17)) to quantify the growth 
rate of disturbances and the onset of instability. 
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(b) 




FIG. 2: (a) Maximum growth rate, Umax; and (b) dominant wavenumber, fcmax for G = 0. 


IV. RESULTS AND DISCUSSION 


In this section we discuss the numerical results and their physical interpretations for dif¬ 
ferent flow parameters. In the absence of viscosity contrast convective instability is featured 
at the miscible interface when a heavier fluid is placed above a lighter fluid. This hydrody¬ 
namic instability is broadly known as density hngering (DF) in the literature. How is this 
convective instability modihed with the viscosity contrast between the underlying fluid? 
Here we investigate the influence of the viscosity contrast on buoyantly unstable miscible 
fluids for both t/ = 0 and U ^ 0, respectively, in MV A and §IVB 


A. Effect of viscosisty contrast in the absence of displacement 

Daniel and Riaz [3] presented hxed interface and moving interface methods to compare 
the natural convections {U = 0) when the dynamic viscosity of the upper fluid is more or 
less than that of the lower fluid. With the help of moving interface method they showed that 
the onset of instability is delayed when the dynamic viscosity increases with depth compared 
to the case of viscosity matched fluids. On the other hand, the instability sets in earlier 
when the dynamic viscosity decreases with depth. This is in contrary to the situation of the 
classical viscous hngering instability in neutrally buoyant huids. 

We revisit the problem of moving interface (cf. [3]) by choosing the characteristic viscosity 
hch = fJ'i, where fii corresponds to the viscosity of the less viscous huid. For example, 
hch = hupp (or /Xch = blow) when /Xupp < /xiow (or /xiow < hupp)- In this case the linear 
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function /(c) and the log-mobility ratio R are defined as, 


/(c) = 


1 C, if /^upp ^ hlo 

C, if /^upp ^ hlo 


i? = In ( — 
l^i 


(18) 


respectively, where Hm corresponds to the viscosity of the more viscous fluid. Linear stability 
analysis is performed for R = 1 and t/ = 0. The temporal evolution of the maximum growth 
rate, cXmax, and the dominant wave number, fcmax, of the perturbation quantities are shown in 
£gure|^ Figure l^a) depicts that the onset happens at the earliest for the classical DF case 
compared to the situations when the viscosity of the two fluids are different. Moreover, the 
temporal evolution of the growth rates obtained for different viscosity contrasts are visually 
indistinguishable. The physical explanation for different onset time of instability in flow 
with or without viscosity contrast in buoyantly unstable miscible fluids can be presented in 
terms of the instantaneous vorticity perturbation field as discussed by Daniel and Riaz [3]. 
From hgure |^b) it is identihed that for all time the most unstable wave numbers are larger 
in the absence of the viscosity contrast as compared to the situation when the viscosity of 
the two fluids are different. Thus we conclude that for [/ = 0 viscous force induces stability 
to the classical DF instability in a vertical porous medium. 


From §IIB| it is clearly observed that all the characteristic scales are derived using the 
characteristic viscosity. Daniel and Riaz |3] chose /ich = /Uupp for both the more and less 
viscous upper fluid, so that the viscosity-concentration relation ([^ takes the form /i(c) = 
gij(i-c)^ where R = ln(/iiow/hupp)- Therefore, R > 0 (or i? < 0) corresponds to the less (or 
more) viscous upper fluid. Such a choice of the characteristic viscosity generates different 
length, time and velocity scales for the respective problem of more or less viscous upper 
fluid. Writing the characteristic length, time, velocity and dynamic viscosity, corresponding 
to R > 0 as and /i)^, and those for R < 0 as Kh” /^ch’ respectively, 

we have 


Fch/F/h = Kh/Ki = LJLX = 1/a, tJtX = 1/ 


a 


(19) 


Thus comparison of the onset of instability and fingering dynamics between these two 
cases should be performed by suitably choosing the characteristic viscosity. In the present 
analysis ni is chosen as the characteristic viscosity irrespective of whether the upper fluid is 
more viscous or less viscous than the lower fluid. Using such a viscosity scaling the length, 
time and velocity scales of the two fluid flow problems corresponding to a more or less viscous 



a 


(i) 

/ich R fJ-ic) Vel. Len. Time 


(ii) 

fich R m(c) Vel. Len. Time 


10 (> 1) //upp 2.3 U L t /Xupp 2.3 U L t 

0.1 (< 1) //low 2.3 U L t //upp -2.3 U/a aL a^t 

TABLE I: The log-mobility ratio (i?), dimensionless velocity (Vel.), dimensionless length 
(Len.) and dimensionless time (Time) corresponding to two different viscosity scales are 
shown for a given set of dimensional valnes. Here U = L = L/ and t = t/r^, 

with " being the dimensional valne of the respective variables. 

ffnid at the top remain the same. For a given set of dimensional valnes of the displacement 
velocity, domain length and time, respective dimensionless valnes corresponding to /ich = 
/^upp mm and /ich = are presented in table It is identified that corresponding to 
l^ch = /Tupp, the dimensionless valnes obtained for a more viscous fluid at the top are different 
from those when the less viscous fluid is at the top. A simple rescaling of the dimensionless 
displacement velocity, length and time of the problem is essential to compare between the 
cases of more {R < 0) or less (R > 0) viscous fluid at the top as shown in table However, 
Daniel and Riaz |3] used U, L and t as the dimensionless values for both R > 0 and R < 0 
(see figures 4, 5 and 12a, etc. of 0). This rescaling can be avoided with /ich = l^i- 

In order to understand the influence of viscosity scaling we consider the problem of 
viscous stabilization of buoyantly unstable miscible layers, i.e. a = /iiow/h-upp < 1, such 
that /ich = /iupp results /i(c) = with R = ln(/iiow//rupp) < 0, while /ich = fiiow 

corresponds to /i(c) = with R = ln(/iupp//iiow) > 0. Therefore, suitably choosing the 
length, time and velocity of the flow problems, one would expect to identify the same 
results for R < 0 and R > 0. In order to illustrate this fact, we choose U = 0 and 
a = 0.1 so that \R\ ~ 2.3 and perform DNS using a Fourier pseudo-spectral method [20] 
to support our theoretical analysis. The obtained numerical results are depicted in figures 
ga-c). Figure g a), which corresponds to the viscosity scaling applied by Daniel and Riaz 
0, i-e. /i(c) = e-2'3(i-'=), shows that the diffusive interface features fingers. On the other 
hand, DNS results corresponding to /ich = /^low (i-e. /i(c) = depict pure diffusive 

expansion of the miscible interface (see figure gb)). We also perform DNS corresponding 
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(a) (b) (c) 



FIG. 3: Spatial distribution of the solute concentration for G = 0: (a) /i(c) = e at 

t = 2 X 10^, (b) /i(c) = at t = 10^ and (c) /i(c) = q,2^ ^ 10^/10^ = 10^. 



FIG. 4: (a) Maximum growth rate, Umax, and (b) dominant wavenumber, fcmax for U = 1. 


to /i(c) = ^ith rescaled length and time scales according to the above-mentioned 

relations (see equation 19). Spatial distribution of the solute concentration at a^t = 10^ 
is shown in £gure|^c), which is identical to figure |^b). From linear theory we identihed 
that the dynamics of the systems for more and less viscous upper fluid are indistinguishable. 
The analogous results can also be shown in the nonlinear regime through DNS. Thus we 
conclude that the comparative study presented by Daniel and Riaz |3j is inappropriate, 
since the length, time and velocity scales associated with the problems related to more or 
less viscous upper fluid are different in their study. 
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B. Influence of the displacement velocity 


Next, we consider the displacement of the lower fluid by the upper one. The influence of 
both the stable and unstable viscosity contrasts are discussed through LSA as well as DNS. 
Figure]^ shows the temporal evolution of cimax and fcmax for the same parameters as those of 
hgurej^ except for 17 = 1. It is identihed that the dynamics of buoyancy induced instability 
in viscosity matched fluids remains unaffected with the dimensionless displacement velocity, 
hence the onset of instability are the same for t/ = 0 (dash-dotted line in figure [^a)) and 
U =1 (dash-dotted line in hgurej^a)) when /r(c) = 1. Figure |^a) depicts that instability 
sets in earlier when the less viscous heavier fluid at the top displaces the more viscous 
lighter fluid at the bottom. The unfavourable viscosity contrast coupled with buoyancy 
force enhances the instability, which is readily evident from the fact that for all time the 
growth rate of the perturbations for /r(c) = e^~'^ are larger than the respective values when 
/r(c) = 1 (see hgurej^a)). On the other hand, for /r(c) = the favourable viscosity contrast 
acts against the instability induced by the buoyancy force and the displacement becomes 
stable (see hgurej^a)). Figure |^b) illustrates that at a given time and for the parameter 
values scanned here, the most unstable wave number fcmax is the largest for /i(c) = 
and smallest for /r(c) = e^. In other words, at a given time /Cmax increases with a. Such 
influences of the viscosity contrast on the stability of buoyantly unstable miscible fluids 
are similar to those in neutrally buoyant fluids, i.e. in the case of VF instability [ 8 ]. In 
summary, for neutrally buoyant as well as buoyantly unstable miscible fluids an unfavourable 
viscosity contrast enhances the instability, whereas a favourable viscosity contrast weakens 
the instability, when the upper fluid displaces the lower one. 


In ( IV A we show, for 17 = 0, the dimensionless length and time should be chosen wisely 
to compare between 7? > 0 and 7? < 0. Here we continue similar analysis for f/ 7 ^ 0. As 
an example we choose U = 0.5, and perform DNS when the dynamic viscosity of the upper 
fluid is more or less than the lower one. Following Daniel and Riaz m we choose fi^h = /^upp 
and compare the dynamics of less viscous upper fluid, i.e. /x(c) = j'ggg figure [^a)) 

with that of the more viscous upper fluid, i.e. /i(c) = j'ggg figure |^b)). Counter¬ 

intuitive results, that the displacement of a less viscous fluid by a more viscous one features 
stronger instability than the displacement of a more viscous fluid by a less viscous one, are 
identihed. Next we take /ich = such that the displacement of more viscous fluid at the 
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FIG. 5: Spatial distribution of the solute concentration: (a) /i(c) = ^\U = 0.5 at 

t = 2 X 10^, (b) /r(c) = e-2'3(i-^), f/ = 0.5 at t = 2 x 10^, (c) /i(c) = G = 0.5 at t = 10^ 
and (d) /r(c) = all = 0.5 x 10 = 5 at a'^t = 10^/10^ = 100. 


bottom by less viscous fluid at the top is represented by /i(c) = and the corresponding 
spatial distribution of the solute concentration is presented in £gure|^c). It depicts that the 
miscible interface features only diffusive expansion, which was also mentioned by Manickam 
and Homsy HI through LSA. Further, we perform numerical simulations for /i(c) = e 
with rescaled dimensionless velocity, length and time as mentioned in table [Tj The result 
obtained from DNS is depicted in £gure|^d), which is indistinguishable from figure [^c). 
This signihes the essence of an appropriate viscosity scaling while comparing the influence 
of more or less viscous fluid at the top on the dynamics of a buoyantly unstable miscible 
interface. 


Discussion of figures |2]j^ illustrates that a simpler and convenient choice of pch is the 
dynamic viscosity of the less viscous fluid, i.e. p./, which preserves the same dimensionless 
length, time, velocity, etc. for more or less viscous upper fluid. In the rest of the paper we 


choose pch = ho such that the viscosity-concentration relation is given by (18). 

The stability scenarios in the phase space spanned by the displacement velocity, U, and 
the mobility ratio, a = piow/hupp, are shown in hgure The parameter space can be 
divided into three distinct regions depending on the stability characteristics. The buoyancy 
and viscosity dominated instability regions are denoted by I and II, respectively, and the 
stable region is represented by region III. For the viscosity matched fluids (i.e. a = 1) 
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FIG. 6: Different stability regions in phase space spanned by the displacement velocity, G, 
and viscosity ratio, a = /iiow/hupp- Regions I and II correspond to the instabilities 
dominated by buoyancy and viscosity, respectively, while the stable region is denoted by 

III. 


the diffusive interface is buoyantly unstable and the growth rate of the perturbations are 
indistinct for all possible values of the dimensionless displacement velocity, U. For a > 1, 
i.e. when a less viscous heavier fluid at the top displaces a more viscous lighter one at the 
bottom, the instability at the diffusive interface is driven by both the viscosity and density 
contrasts, and the instability increases with a as well as U. On the other hand, for a < 1, 
i.e. for viscous stabilisation of buoyantly unstable diffusive interface, the instability becomes 
weak as a decreases or U increases, and hnally becomes stable when U (or a) is larger (or 
smaller) than a certain critical value. 

The influence of the viscosity contrast on the onset of instability is depicted in hgure 
This hgure shows, in the absence of displacement {U = 0), the onset of instability delays with 
the viscosity contrast, irrespectively whether the heavier huid is less viscous or more viscous. 
The present results are consistent with Daniel and Riaz [3] when a more viscous heavier 
huid is placed over a less viscous lighter huid. These authors used /Xupp as the characteristic 
viscosity and showed, for t/ = 0 the onset time increases with mobility ratio a = /xiow/hupp- 
More surprisingly, they found an early onset of instability with stable viscosity contrasts 
(a < 1) and U > 0, compared to buoyantly unstable miscible interface between viscosity 
matched huids, i.e. a = 1 (see hgure 12a of Daniel and Riaz |3]). This result is contradictory 
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FIG. 7: Effect of viscosity contrast on the onset time, tc, for different dimensionless 
displacement velocity, U = 1 (solid line), U = 0.5 (dashed line) and U = 0 (dash-dotted 
line). Sqnare (□) corresponds to /i(c) = and asterisk (*) corresponds to /r(c) = 


(a) (b) (c) 



FIG. 8: Spatial distribntion of the concentration, c, for U = 1: (a) /i(c) = at t = 10000, 
(b) p(c) = 1 at t = 2400 and (c) p(c) = at t = 1000. 

to the resnlt of Manickam and Homsy na, who showed that a buoyantly unstable miscible 
interface can be stabilised by suitably choosing the viscosity contrast and the displacement 
velocity. The present LSA results successfully captures this phenomenon (see £gure|^. The 
onset of instability as a function of the log-mobility ratio, both for more and less viscous 
upper fluid, is depicted in figure It is identified that the instability sets in earlier with 
increasing R when a less viscous fluid displaces a more viscous fluid. On the other hand, 
during the displacement of less viscous fluid by a more viscous one, onset time increase 
with R and the displacement becomes completely stable after a threshold value of R, which 
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depends on the injection velocity U. In order to confirm the present linear stability results, 
DNS are performed for the parameter values used in figure]^ The spatial distribution of the 
concentration for displacements with and without viscosity contrast are shown in figure 
This figure depicts that a buoyantly unstable diffusive interface of viscosity matched fluids 
(middle panel in figure becomes stable when the dynamic viscosity of the upper fluid is 
more (left panel of figure |^. On the other hand, the instability becomes stronger when 
a less viscous fluid at the top displaces a more viscous fluid at the bottom (right panel of 
figure . 


C. Non-orthogonality of eigenmodes and transient growth 


Here, we briefly discuss about the possible transient growth of the disturbances to the 
unsteady base state Ch{x,t) = 0.5 erfc(a;/2\/t). Manickam and Homsy [H] presented an 
LSA based on the modal analysis under the assumption of quasi-steady state approximation 
(QSSA). However, the QSSA modes are non-orthogonal P, [TH] and thus these eigenvalues 
do not reveal the exotic transient behavior (Schmid2007). Further, Trefethen et al. [22] 
reported that the transient growth in a stable linearised system has implications for the 
behaviour of the associated nonlinear system. 

To quantify the degree of non-orthogonality of the eigenmodes obtained from modal 
analysis we study the linear stability of the unsteady base state Ch{x, t) with respect to small 
wavelike perturbations of the form 


(c',m') {x,y,t) = {4>c,4>u) {x,t)exp{iky), 


( 20 ) 


where i = a/— 1, k is the non-dimensional wave number in the y direction, and (/)c(a;, f), cf)u{x, t) 
are time dependent concentration and velocity perturbations, respectively. Following the 
standard procedure n E], the linear stability equations in a similarity transformation 
(^, t)-domain can be written as an initial value problem (IVP) for 0c, 


d(t)c (I d'^ , 2V ^ ^ \ , 1 dcfe 


d 


d 


— + Rficf,)— - kH (t>u = kH UR fie) - 


P'{c,) 


( 21 ) 

( 22 ) 


j V ' ' po 

where ^ = x/\/t is the similarity variable. Finite difference approximation of the linearized 


operators followed by elimination of 0„ from equations (21) and (22) yields a nonautonomous 
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system of ordinary differential equations, 


^ = A{k, ti) = 


oo < ^ < oo, 


(23) 


where ti corresponds to the initial time, when perturbations are introduced and A{k, t) is 
the time dependent matrix. 

For the given matrix A{k,t), we would like to have some effective way to determine 
whether one should be concerned about the effects of non-normality. The simplest quantita¬ 
tive approach often used for characterising normality is, k{V) = || 1^||2 I I 2 , the condition 
number of the eigenvector matrix V associated with A{k,t) [ 6 j. Here 11-112 corresponds to 
the standard Euclidean norm. It can be shown that for a normal matrix A{k, t) the condi¬ 
tion number niV) is 1. In order to quantify the potential transient growth of disturbances 


and non-normality of A{kA) associated to the IVP (23), first we compute k{V) and then 
the numerical abscissa and the spectral abscissa, denoted by, a{A) and ri{A), respectively, 
and are dehned as. 


a{A) = max{3fJ(A(^)}, 
rf^A) = max{A(Al -|- Ai^)/2)}. 


(24) 

(25) 


Here A = AikA), A(-) represents the eigenvalue of the respective matrices, 3fJ(-) denotes 
the real part and AA- denotes the transpose of the matrix A. The numerical abscissa ri{A) 
measures the maximum possible instantaneous growth rate corresponding to any initial 
condition as t —)■ 0 [ 21 ]. It is important to note that for a normal matrix a{A) = f]{A). The 
scalar measures A(-) and k{-) of non-normality of the matrix A{k, t) are computed using the 
MATLAB routines eig and cond, respectively. 

The effect of non-normality in terms of the condition number, numerical abscissa and 
spectral abscissa for various flow conditions are summarized in tables [IT] and ES The time 
dependent matrix A is frozen at different time and the computed a{A),ri{A) and k(H) are 
tabulated in table for classical VF of neutrally buoyant fluids and DF of viscosity matched 
fluids for a given wave number A; = 0.1. For the case of VF the displacement velocity is taken 
as the characteristic velocity, such that U = 1, and the log-mobility ratio is R = 1. Next, we 
consider the influence of viscosity contrast on buoyantly unstable miscible fluids for U = 0 


as well as [/ 7 ^ 0. Table HI compares between two cases: (i) /i(c) = e^^p{c) = c,U = 0 and 
(ii) /i(c) = e^“'^,p(c) = c, f/ = 1. Since k(V) is very large for all the cases discussed here. 
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(i) (ii) 

U = 0,/x(c) = l,/ 5 (c) = c U = l,/u(c) = e^“'^,/ 9 (c) = 1 

(DF in viscosity matched fluids) (VF in density matched fluids) 


to 

a{A) ri{A) 

k{V) 

a{A) v{A) 

k{V) 

0.1 

-4.9626 -2.4719 

3.2121e-h24 

-4.9658 -2.4708 

3.9351e-h24 

0.5 

-0.9652 -0.4730 

7.5720e-F24 

-0.9677 -0.4724 

9.2730e-h24 

1 

-0.4671 -0.2243 

4.7560e-F25 

-0.4692 -0.2239 

4.6768e-h24 

5 

-0.0738 -0.0291 

4.8664e-F24 

-0.0750 -0.0290 

8.9276e-h24 

10 

-0.0276 -0.0069 

2.6092e-F25 

-0.0285 -0.0070 

6.4083e-h24 

20 

-0.0069 0.0023 

1.8240e-F25 

-0.0075 0.0022 

9.5316e-h24 

30 

-0.0013 0.0043 

3.6163e-F25 

-0.0017 0.0043 

2.5617e-h25 

50 

0.0020 0.0050 

4.4731e-h24 

0.0018 0.0050 

5.2046e-h25 


TABLE II; For a given wave number k = 0.1, the variation of spectral abscissa C({A), 
numerical abscissa ri{A), and the condition number of eigenvector matrix k{V), at different 
frozen time to: (i) DF in viscosity matched fluids, (ii) VF in density matched fluids. 


a substantial non-modal growth of the disturbances at early time can be anticipated, and 
this is confirmed from the difference between a{A) and ri{A) during initial period. Table 
[II] depicts that the order of non-normality is almost equal for these two cases. From table 


III it is identihed that for U = 0 both a{A) and r][A) are negative, which signihes that in 


the presence of the viscosity contrast onset of instability is delayed. On the other hand, for 
t/ 7 ^ 0 instability is enhanced. These results are consistent with our observation from LSA 


as well as DNS discussed in OV 


Further, in order to understand the importance of appropriate viscosity scaling on the 
non-normal growth of the perturbations a{A) and t]{A) are plotted in figure]^ Figure [^a) 
depicts that for a given wave number k = 0.05 and D = 0 , a{A) and r]{A) corresponding 
to fi{c) = are identical to those corresponding to /i(c) = The influence of the 

viscosity contrast in the presence of fluid displacement (i.e. 1 / 7 ^ 0 ) on non-normality of 
the linearised matrix A is presented in hgure |^b) for U = 0.5, i? = 0.25 at to = 30. This 
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(i) (ii) 

[/ = 0, /i(c) = e^, p{c) = c U = 1, p{c) = p{c) = c 

(More viscous fluid at the top) (Less viscous fluid at the top) 


to 

a{A) rj{A) 

k{V) 

a{A) 

r]{A) 

<V) 

0.1 

-4.9841 -2.4930 

5.4772e-F24 

-4.9399 

-2.4427 

7.4439e-h24 

0.5 

-0.9851 -0.4927 

5.4482e-F24 

-0.9429 

-0.4457 

1.2781e-F25 

1 

-0.4859 -0.2429 

1.6465e-F25 

-0.4452 

-0.1984 

9.9258e-h24 

5 

-0.0890 -0.0441 

2.4518e-F25 

-0.0540 

-0.0074 

3.1297e-h24 

10 

-0.0450 -0.0200 

2.0467e-F24 

-0.0094 

0.0125 

3.3447e-h25 

20 

-0.0183 -0.0088 

6.1171e-F24 

0.0093 

0.0193 

8.1408e-F24 

30 

-0.0114 -0.0055 

5.5201e-F24 

0.0138 

0.0199 

1.0704e-F25 

50 

-0.0066 -0.0034 

5.3213e-F24 

0.0154 

0.0187 

2.2918e-h25 


TABLE III; For a given wave number k = 0.1, the variation of spectral abscissa a{A), 
numerical abscissa ri{A), and the condition number of eigenvector matrix k{V), at different 
frozen time to: (i) the influence of a stable viscosity contrast in the absence of 
displacement, (ii) the influence of an unstable viscosity contrast in the presence of 

displacement. 


hgure depicts that the numerical (spectral) abscissa corresponding to /i(c) = is larger 

than that corresponding to p,{c) = Thus, we conclude that by choosing an appropriate 
characteristic viscosity one can lead to the same non-normal growth of the perturbations 
associated to the problem of a buoyantly unstable miscible interface both with the stable 
and unstable viscosity contrast when U = 0. On the other hand, in the presence of fluid 
displacement the instability is stronger when the less viscous fluid at the top displaces the 
more viscous fluid at the bottom. These observations are consistent with the LSA as well as 
the DNS presented in §IV| and could be captured only through the scaling analysis discussed 
in the present paper. 

To summarize, it is observed that the non-normality of the linearised matrix in the 
study of hydrodynamic instability driven by buoyancy or viscosity or both is of signihcant 
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FIG. 9: Numerical abscissa (Q) and spectral abscissa (□): (a) 
k = 0.05, t/ = 0, = 1, to = 50, (b) k = 0.05, U = 0.5, R = 0.25, R = 30. 

Viscosity-concentration relation is colour coded; /i(c) = (red) and /i(c) = 

(black). 


importance at early time. Although the frozen time approach to measuring the degree of 
non-normality of the time dependent matrix A through r]{A) and/or k{V) only provides a 
crude approximation, it can be handy to obtain an insight about possible non-modal growth 
of the disturbances. The transient growth of perturbations in a non-autonomous system 
can be determined efficiently through the propagator or matricant approach [91 [H] or the 
direct adjoint looping (DAL) analysis [HE], which is beyond the scope of the current study. 
To determine the optimal perturbation leading to the instability in such cases is the topic 
of ongoing research and it is strongly believed that the importance of the scaling analysis 
discussed in this paper can also be observed in the optimal perturbations. 


V. CONCLUSION 

We numerically investigate the influence of viscosity contrast on buoyantly unstable mis¬ 
cible interface in vertical porous media using an LSA as well as DNS. We show that in the 
absence of displacement a buoyantly unstable viscous miscible interface is the least stable 
when the viscosity of two fluids are equal, compared to the variable viscosity interface. In 
this case instability sets in at the same time for both less and more viscous upper fluid. On 
the other hand, a less viscous heavier fluid displacing a more viscous lighter fluid features an 
earlier onset than when the more viscous heavier fluid displaces the less viscous lighter fluid. 
We also show how a suitable rescaling of the dimensionless length, time and the displace¬ 
ment velocity can reproduce the results of Daniel and Riaz [3] from the present analysis and 
vice-versa. Thus the importance of an appropriate scaling analysis in fluid mechanics prob- 
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lems is presented by investigating the influence of viscosity contrast on buoyantly unstable 
miscible fluids in vertical porous media. The principal aim of an LSA is to obtain the onset 
of instability accurately and to predict the optimal perturbation that leads to the instability. 
Non-modal analysis to determine the optimal growth in buoyantly unstable miscible fluids 
with viscosity contrast has been undertaken for further study. 
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